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ABSTRACT 

We present a numerical code intended for calculating stellar evolution in close 
binary systems. In doing so, we consider that mass transfer episodes occur when the 
stellar size overflows the corresponding Roche lobe. In such situation we equate the 
radius of the star with the equivalent radius of the Roche lobe. This equation is handled 
implicitly together with those corresponding to the whole structure of the star. We 
describe in detail the necessary modifications to the standard Henyey technique for 
treating the mass loss rate implicitly together with thin outer layers integrations. 

We have applied this code to the calculation of the formation of low mass, helium 
white dwarfs in low mass close binary systems. We found that the global numerical 
convergence properties are fairly good. In particular, the onset and end of mass transfer 
episodes is computed automatically. 
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1 INTRODUCTION 

Stellar evolution in binary systems is a classical topic of 
stellar astrophysics. The study of stellar evolution in this 
kind of systems has great importance, since an important 
fraction of know stars belong to binary systems. Inasmuch 
as, if we want to build a stellar evolution theory able to 
account for the observations, it is essential to study stellar 
evolution in binary systems. 

According to the standard model, the initially most 
massive star, usually named the primary evolves faster and, 
as consequence of its nuclear evolution, begin to innate. 
If stars are sufficiently close each other (close binary sys- 
tems, hereafter CBS) the primary will overflow of the "Roche 
lobe" . This lobe is found by considering the problem of a bi- 
nary system in which the orbits are circular, and rotation of 
both objects is synchronized with orbital motion. Therefore 
we have the Largange solution to the restricted three body 
problem, considering that each gas particle is of infinitesi- 
mal mass. The Roche lobe corresponds to an equipotential 
surface that surrounds the stars and represents the limiting 
volume that can be attained by the stars of the pair before 
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the beginning of mass transfer either, onto the companion, 
and/or outside the system. From this moment on, the evo- 
lution of each star is markedly different from the one they 
would have in isolation. 

According to Kippenhahn & Weigert (1967) mass trans- 
fer may begin during core hydrogen burning (case A) and 
also during central contraction after central hydrogen ex- 
haustion (case B). If it begins after core helium burning it 
is commonly referred to as case C (Lauterborn 1970). Mass 
transfer can occur in conservative conditions (total mass and 
orbital angular momentum conserved) or with mass and/or 
angular momentum losses from the system. These losses 
can drastically affect the evolution of the star in a CBS. 
For reviews on stellar evolution in binary systems, see, e.g., 
Paczyhski (1971), Iben (1985), Iben & Livio (1993), and 
Taam & Sandquist(2000). 

In view of the large number of situations that can be 
encountered, it is not surprising that CBSs are considered 
as favorable places for the occurrence of a number of phe- 
nomenologies with high astrophysical interest: X-ray sources 
(van der Klis 2002) , low mass helium white dwarfs (Kippen- 
hahn, Kohl & Weigert 1967), planetary nebulae with binary 
nucleus (Iben & Livio 1993), binary radio pulsars (Bhat- 
tacharya & van den Heuvel 1991), supernova progenitors 
(Podsiadlowski, Joss & Hsu 1992). In this last case, binary 
evolution naturally accounts for the evolution of the progen- 
itor of SN 1987A. This object was a blue supergiant, a fact 
difficult to explain in the frame of single stellar evolution, 
but not in binary evolution (Podsiadlowski, et al. 1992). 
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In view of the importance of this line of research, in our 
Observatory we have begun the study of CBSs. To this end, 
we have developed a code to compute stellar evolution in 
CBSs, including the self consistent calculation of the mass 
transfer rate M (hereafter MTR). The main purpose of the 
present paper is to describe the numerical techniques em- 
ployed in our program for the computation the evolution 
of stars in a CBS. This code is based on a modification of 
the Henyey technique presented in Kippenhahn, Weigert & 
Hofmeister (1967) to solve the set of difference equations of 
stellar evolution with the inclusion of hydrodynamical ef- 
fects and the MTR. As a first application of this new code 
to the case of CBSs that lead to the formation of helium 
white dwarfs. 

As stated above, the main characteristic of our code 
is that, during mass transfer episodes, it finds the MTR in 
a fully implicit way by means of a single iterative proce- 
dure. This is in contrast with the usual procedure employed 
by other researchers (see, e.g., Podsiadlowski, Rappaport, 
& Pfahl 2001; Wellstein, Langer & Braun 2001). In these 
works, in the construction of a stellar model a double it- 
erative procedure is applied: a MTR is estimated, and the 
stellar structure is relaxed for this MTR value. Then, with 
this new structure, a new MTR is computed and the stellar 
structure is relaxed again. This is repeated until consistency 
is achieved. In calculating the MTR it is usual to employ the 
physical treatment presented by Ritter (1988) which gives 



M = —M exp 



Rl-R 
Hp 



(1) 



where Mo is the MTR for a star that exactly fills the Roche 
lobe, R L is the equivalent radius of the Roche lobe (see 
below), R is the stellar radius and Hp is the photospheric 
pressure scale height. Notice that the double iterative loop 
is necessary to get numerical stability because of the steep 
dependence of the MTR upon the sizes of the star and the 
Roche lobe. 

As the main goal of this paper is to present the numer- 
ical scheme we have devised and not to produce start - of - 
the - art evolutionary models here, for simplicity, we shall 
simply impose Rl = R- In doing so, we shall get a good 
global computation (but not the details of the beginning 
and end) of each mass transfer episode. 

The rest of the paper is organized as follows: In Sec- 
tion 2 we briefly describe the differential equations of single 
stellar evolution we have to solve. In the following Section 
3, we describe in detail the Henyey technique and also the 
code employed to compute binary stellar evolution. We de- 
scribe our treatment for the orbit of the binary, including 
the modification we have considered for total mass and an- 
gular momentum losses from the system in §4. In Section 5 
we discuss the critical problem of how to manage the outer 
boundary conditions that, in turn, are essential in determin- 
ing the MTR during mass transfer episodes. Then, in Section 
6 we discuss some numerical techniques we have applied, es- 
pecially the rezoning of the stellar models. In Section 7, we 
present the results for three different binary systems that 
lead to the formation of helium white dwarfs. Finally, in 
Section 8 we comment on our results, further improvements 
we plan to incorporate to the code and the astrophysical 
problems we plan to study with. 



2 THE EQUATIONS OF STELLAR 
EVOLUTION 

Here, we briefly summarize the equations of stellar evolu- 
tion to be solved. As usual, we consider spherically symmet- 
ric objects, neglecting rotation and magnetic fields. In these 
conditions, the equations of stellar structure are (for deriva- 
tion of these equations see, e.g., Clayton 1968, Kippenhahn 
& Weigert 1990. For a detailed treatment of hydrodynamic 
stellar codes, see, e.g., Kutter & Sparks 1972): 



i) the Euler equation of fluid motion 

dv ^ dv _ 1 dP G m r 
dt dr p dr r 2 

ii) the definition of velocity 
dr 

di= V > 

iii) the equation of mass conservation 

drrir . 2 
— = 47rr p, 

iv) the equation of energy balance 



Olr . 2 / 

— = 4nr p I e 



■£„ -T 



dS\ 
dt) 



(2) 



(3) 



(4) 



(5) 



v) the equation of energy transport for the radiative case 
dT 3-7T up l r 

~dr = ~4acT3 47r 2 r 2 ' ( ' 

and 

vi) the equation of energy transport for the convective case 

or rap 

Or ~ conv P dr ■ ( ' 

We employ the Schwarzschild criterium for the onset of con- 
vection. The symbols have their standard meaning. 



3 THE HYDRODYNAMIC CODE FOR 
BINARY STELLAR EVOLUTION 

In order to incorporate the specific phenomena occurring in 
binary evolution we have to make supplementary assump- 
tions apart from those quoted in the previous section. We 
shall handle the members of our system as spherical ob- 
jects, neglecting the departure from spherical symmetry of 
the equipotentials (e.g., the Roche lobe) and its evolution- 
ary consequences. Moreover, we shall consider that the ob- 
jects move along circular orbits, and that they influence each 
other only through gravitational attraction (we neglect irra- 
diation). 

As usual we shall consider the problem in Lagrangian 
coordinates, considering £ the independent variable, defined 

as 



C = ln 



(8) 



Radius, pressure, and temperature are handled by means of 
logarithmic transformations 

p = \nP, 
9 = InT, 
x — lnr, 
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whereas l r , v are considered linearly. 

For simplicity, we have written the difference equations 
in a centered fashion. It means that we have chosen to write 
a generic differential equation 



(9) 



J ' J - F ( x j+i/2,yi,j+i/2, ■ ■ ■ ,115,3+1/2) 



(10) 



■ 2/sj; 

as a difference equation 

Vi,3 + 1 
Xj-\-\ — Xj 

where fJj+i/2 = (%+i + %)/2, being rj any quantity. The 
second subindex j indicates the shell of the star for which 
the difference equation is written. The results obtained in 
this way are completely satisfactory for our present pur- 
poses. However, this would have not been the case if we 
were studying stellar objects suffering shock waves some- 
where. For calculations of more violent phases we shall have 
to include non - centered equations together with some dis- 
sipative effect (e.g., artificial viscosity, see Kutter & Sparks 
1972). Temporal derivatives have been written in the stan- 
dard backwards differenced form. 

In problems with variable stellar mass we need to be 
careful at computing temporal derivatives at constant mass. 
We have found it very convenient to rewrite the derivative 
operator as 



(11) 



In this formulation we get the dependence of these deriva- 
tives with the MTR. This is important because we shall treat 
the MTR as a new variable to be relaxed in our code. Then 
we have to design a Henyey scheme for relaxing the internal 
structure of the star, the total luminosity, the effective tem- 
perature and (in the case of the occurrence of mass transfer) 
the MTR. 

As usual, the calculation of the corrections is based in 
the solution of the lineal system 



d 


d 






d 


dt 


~ dt 









H- 5 = h 



(12) 



where H is a matrix with a structure similar to those cor- 
responding to the cases of single star evolution or codes in 
which the MTR is not treated self consistently. It has blocks 
on the diagonal, but now, an important difference appears: 
now there is a non vanishing column due to the derivatives 
of the structure equations with respect to the MTR. The 
vector 8 represents the corrections vector (see below) and h 
is the inhomogeneity vector that must vanish for a correct 
stellar model. 

Here, some words are in order about the necessity of 
such a deep modification of the standard numerical treat- 
ments of binary evolution. Savonije (1978) has treated MTR 
as another variable to be relaxed by iterations but only con- 
sidered it to be relevant in the outer layers integration where 
he included most of the change in luminosity due to mass 
transfer (See Section 5 for details). This is so only in the case 
of excluding a very thick portion of the star from the grid. In 
this way we are neglecting the temporal derivatives, and so, 
the inertia of the outermost layers of the star. Other scheme 
intended for the computation of single and binary hydro- 
static stellar evolution has been presented by Ziolkowski 
(1970), who incorporated an iterative procedure for getting 



the MTR. However, the outer envelope integrations are per- 
formed in a very thick portion of the star, and because of 
the reasons quoted above, this is not entirely adequate for 
our purposes 

In order to perform a realistic simulation of the evo- 
lution of the stars, it is quite desirable to handle most of 
the mass of the star inside finite - difference part of the 
Henyey scheme. Notice that there are many relevant situa- 
tions in which the outer layers structure is a key ingredient 
in determining the evolution of the stars. Regarding binary 
evolution, we may cite for example that in order to compute 
the mass transfer episodes occurring in low mass pre - white 
dwarf stars we do need to consider very thin outer layers 
integrations (see below for more details). 

Consequently, we have to consider an outer layers inte- 
gration in a portion of the star so thin that a fraction, or 
even most, of the luminosity change (see Section 5 for de- 
tails) due to mass loss during mass transfer episodes occurs 
inside the finite - differenced portion of the model. Thus, 
the form of the equation of energy conservation we have to 
consider is 



dl r 



-e ~t( — 

nuc £ " \dt 

where, from Eq. (8), we have 



+ 



dS_ 



(13) 



dt 



M 
M 



(e-«-l). 



(14) 



Thus, we are forced to consider M as a extra Henyey un- 
known. 

Notice that, for consistency, we should consider the to- 
tal mass value of the star as another unknown to be found 
during the relaxation process by means of the equation 



M = M prev + M At, 



(15) 



where M prev is the mass of the previous model. Because of 
the transformation of variables we have considered, the to- 
tal mass value occurs in most of the equations of structure. 
Thus, the matrix we shall have to handle in finding the cor- 
rections will have the usual blocks on the diagonal but also a 
non - vanishing column. We shall describe in the rest of this 
section the technique we employed to solve the numerical 
problem. 

The first block of the matrix includes the four boundary 
conditions equations that link the values of the structural 
quantities at the first meshpoint with the effective tempera- 
ture, luminosity and MTR value. These are constructed by 
meand of linear interpolations among the outer layers inte- 
grations. The resulting block has the form 



dhi 



dhi dh± 

8M ~8~L 

8h 2 8h 2 

dM 8L 9T eff 



8T e 

an 



ff 



dhi 




dhi 


dhi ' 


dxx 


dh 


dpi 


801 


8h 2 


dh 2 


dh 2 


dh 2 


dxi 


dh 


dpi 


801 



8hp 
8M 



8h 5 
8L 



8h 5 



8T C 



ff 



8h 5 
dxi 



8h 5 
dh 



8h 5 
dpi 



<>h r , 
80i 





8M 




SL 




ST eff 




Sxi 




Sh 




8pi 




80! 
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-hi 
-hi 



-h 5 



(16) 



where hi , i = 1, • • • , 5 are functions known from linealizacion 
of stellar structure equations. Hereafter in this section, L is 
log(L/LQ) and T e // the logarithm of effective temperature. 
We define vectors Ui, Vi, Wi, Xi, i = 1, 



5M 
SL 

STeff 

5h 



Ui Vi Wi Xi 

u 2 v 2 w 2 x 2 



U 5 V 5 W 5 X 5 J 



5vi 
Spi 
S0i 
1 



(17) 



For the first block are U% = 0, i = 1, • • • , 5 and 5vi = 0. 
From Eq. (16) and Eq. (17) we can write 



dhi dhi 

SM dL 

dhf dh 2 

dM dL 



dhi 



dhi 
dx\ 

Sh2 

dT eff dxi 



9T e 

an 



ff 



dhp 
dM 



dh 5 
dL 



dh 5 



eft 



9fei_ 

dpi 
dh 2 
dpi 



9fe 5 
dpi 



dhi 
' 98i 
dh 2 



Ohr, 

' d0i 



dh 5 
dxi 

-hi 
-hi 



9hi 
9li 
9h 2 
9h 



dh 



Ui Vi Wi X! 

U 2 V 2 Wi x 2 



U 5 V B W 5 X 5 



(18) 



and from this we can obtain vector components Vi, Wi, Xi, 
i = 1, • • • ,5. In a similar way we can generalize this mecha- 
nism for any interior block. If we propose 

5M = A n 5v n + B n 5p n + C n 50 n + D n , (19) 

we find, generalizing the equations of the interiors blocks, 



A n — A n -lU5n-4 + B n -\U$ n -Z + C n -\U<bn-1 

B n = A n -lV5n-4 + B n -\Vsn-3 + C n -1 Vsn-2 

C n = A n -lW5n-4 + B n -lWs n -3 + C„_l Wr,n- 2 

D n — A n -\Xr m -4, + B n -lX5n-3 + Cn-lX^n-1 + D n -1- 



(20) 



The equation that leads to intermediates blocks are in the 
form 

(JO' (JO' 
(SVnOti + 8p n (3i + 89„^i) + 5x n +l ^ 5l n +l 



-5v„ 



+1 



dgi 

dv n+ 



5p n -i 



dgi 
dp n +i 



56 n 



dx n+ i 
dgt 



+i 



dln+l 

£i, (21) 



i = l,2,"-,5 

n = 2,3,---,N -2. 



The gi functions are know from linealizacion of stellar struc- 
ture equations in the intermediate shells, and ai,/3i,fi, f; are 
defined as 



ft — B 22t- + Vr 1 _)_ T/ K 9£i 4. 9gi_ 



Sin ^ 9p„ 



7< = Cn fg + W 5 n-1 



Written in matrix form Eq. (21), with 



5v n 

Spn 
56n 
5Xn+l 
5ln + l 

5Vn+l 
5pn+l 
56n+l 
1 



U'm + l Vsn + l Wsn + l Xr,n + 1 
U$n+ 2 Vsn+1 W$n+1 X 5n + 2 



Usn+5 Vsn+5 Wsn+5 X$n+5 



(22) 



and substituting Eq. (22) in matrix form of Eq. (21), we 
obtain an equation that allow us to find the components for 
vectors U, V, W,X: 



ai Pi 71 
on Pi 72 



dgi 



9x 



ra+l "'n + 1 



9g5 



«« to 75 etr ± a^T ^ 

t^Sn+l Vs„+1 Wsn + l X 5 n + 1 
Ur,n+ 2 Vs n +1 W;,n + 2 X^n+2 



Ur,n + 5 V5n + 5 W5n+5 ^5n+5 . 

-6 

<2 



9gi 


9gi 


dgi 


9v n+ i 
9g 2 


9p„+i 
9g 2 


9«„ + i 
992 


9v n+ i 


9p„+i 


9e„ + i 



dv„ 



dp„ 



dgs 



-6 



(23) 



In the last block we have 8l n = 5x n = (5?j n = 0, thus, 
rewriting the previous matrix equation in the usual way and 
proposing, like we have done in this derivation 



5v N -i 




5pN-l 









V5N-4 W5N-4 X5N-4 

VbaT-3 X5N-3 
V5JV-2 WsN-2 Xr,N-2 





SpN 




56 N 




1 



(24) 



we can compute all components of vectors U,V,W,X, and 
then all the corrections. Finally, at the surface we have to- 
tal luminosity, effective temperature and the MTR for each 
model. 

As stated above, in this work we assume that mass 
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transfer occurs if stellar radius R overflows the effective ra- 
dius Rl of the Roche lobe, calculated with the formula of 
Eggleton (1983), 



Rl 
a 



0.49- 



..2/3 



0.6 ■ g 2 / 3 + In (1 + q 1 / 3 



(25) 



where a is the orbital separation and q = M/M2 the mass 
ratio of the binary components. 

When a star losses mass, either by a stellar wind in the 
case of isolated star or by transfer in CBS evolution, not 
only the external layers, but all stellar structure is affected 
by this phenomenon. An outer layer that was internal, af- 
ter a certain time interval, it will be in the stellar surface. 
Therefore, we need to consider in the code this important 
effect. We have found it enough to simulate this outflow by 
means of a linear interpolation in the abundances per mass 
fractions. Then we should check the meshpoint distribution 
for consistency. 

It is important to stress here that in doing so, we are 
introducing noise in the chemical profile of the donor star. 
Presumably, this is responsible (at least in part) for the noise 
in the MTR vs. time relations shown below for the specific 
systems we have computed. This is especially noticeable in 
the case of donor stars with convective envelopes. Evidently, 
a better strategy should be to employ a moving grid tailored 
to follow the outwards motion of the stellar layers. We plan 
to incorporate this refinement in the near future. 



4 THE ORBITAL EVOLUTION 

As a first approximation, we can perform conservative bi- 
nary evolution calculations. In this case we consider total 
mass and orbital angular moment as constants. However, 
we expect the occurrence of total mass and angular mo- 
mentum losses to be very important in some astrophysical 
situations of interest. Thus, we have included in our code, 
these phenomena. 

From the definition of total angular orbital moment, 
and using Kepler's third law, we can write 



2j + A/ 



\M + M 2 M J \M + 



Ah 



2 
W 2 



)(26) 



where J/ J represents energetic losses due to different pro- 
cesses, M and M2 the masses, and M, M2 the MTRs of the 
lossing - accreting stars. 

If we consider non - conservative mass transfer, this 
causes an angular loss from the system. We follow the 
method developed for Rappaport, Joss & Webbink (1982) 
and Rappaport, Joss & Verbunt (1983). It is specified by 
two free parameters, the fraction (3, of the mass lost by the 
primary star that is accreted by the secondary star 2 and 
the specific angular momentum, a, of matter lost from the 
system in units of 2na 2 /P, so that 



6J = aSM(l- /3) 



(27) 



where SM is an incremental mass lost by the donor star, 8 J 
the incremental angular momentum of the matter lost by 
the system, and a and P are the orbital semi major axis and 
period, respectively. We assume that the orbit is always well 
approximated by a circle of radius a (where a is a function 



of time). If we consider only angular momentum losses due 
to the ejection of matter from the system, the last equation 
can be rewritten as the differential equation 



5J = q(1 - j3)y/G (M + M 2 ) a 8M. 



(28) 



Using Kepler's third law to eliminate P, and combined 
with an expression for the total systemic angular momen- 
tum, 



J = 



Ga 



M + M 2 



-MMo 



(29) 



Eq. (28) constitutes a differential equation for J (or a) as 
a function of M (we have neglected the rotational angular 
momentum of the components stars). Then, we write the 
loss of angular momentum for the matter ejection as 



dlnJuE _ «(1 - P)yJ(M + M 2 )Ga 



M. 



(30) 



dt J 
Angular momentum loss due gravitational radiation is calcu- 
lated according to the standard formula (Landau & Lifshitz 
1959) 



d In Jgr 
dt 



32GV Ml t 
5c 5 a 4 



(31) 



where M tot = M + M 2 and fx = MM 2 /(M + M 2 ), G and 
c are the gravitational constant and vacuum speed of light, 
respectively. 

To calculate the angular momentum loss due to mag- 
netic braking, we use the prescription of Rappaport, Joss & 
Verbunt (1983), which is based on the magnetic - braking 
law of Verbunt & Zwaan (1981), 



= -3.8 x IO- 30 MRWdyn cm. 
dt y 



(32) 



where ui is the angular rotation frequency of the donor star, 
assumed to be synchronized whit the orbit. We include full 
magnetic braking when the star has a sizeable convective 
envelope embracing a mass fraction > 0.02. If we substitute 
Eqs. (30) - (32) in (26), and consider in this last equation 



M 2 = -/3M, 



(33) 



in view of the definition of (3, we obtain a differential equa- 
tion for the orbital separation, which has no analytical so- 
lution. 

Let us remark that the semiaxis of the orbit is depen- 
dent upon the MTR and so, the equivalent radius of the 
Roche lobe, which in turn, has been assumed to be the ra- 
dius of the star. Thus, during mass transfer episodes, to be 
consistent with the iterated value of MTR, we need to per- 
form an orbital integration for each iteration. As the equa- 
tions for the orbital evolution are well behaved, we decided 
to handle them with a standard Runge - Kutta technique 
(Press, et al. 1992). This is very fast and represents a tiny 
increment of the total numerical effort. 



5 HANDLING OUTER LAYERS 
INTEGRATIONS 

Here we shall treat the problem of handling the outer layers 
integrations adequately for the case in which there is mass 
transfer. If not, outer layers are managed as in Kippenhahn, 
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et al. (1967). We refer the reader to that paper for further 
details in this simple situation and only describe here its 
modifications to account for the occurrence of mass transfer. 

In the outer layers, we integrate the equations of struc- 
ture neglecting temporal derivatives and the velocity but 
taking the acceleration of these layers into account. Then, 
we have a system of ordinary differential equations. In par- 
ticular, the equation for the luminosity we consider is (see 

Eq 13)* 



dl r = _ r d? 
dm r dt 



dS_ 



(34) 



In the case of constant mass evolution, it simply implies 
that the luminosity is constant in the outer layers. However, 
in the case of mass transfer/loss episodes, it modifies the 
profile of luminosity in a sizeable way. For example, for ra- 
diative envelopes, which have an entropy increase outwards, 
this equation predicts a drop in the luminosity. This is an 
important effect to be considered (see the previous sections). 
If we take a little amount of mass Moli in the outer lay- 
ers integrations (e.g. Moli/M < 1CT 6 ), it is present in the 
outer layers integrations as well as in the finite - differenced 
portion of the star. 

In order to perform an adequate treatment of the outer 
layers of the mass losing component of the pair we need to 
be careful. Because of the way we have planned our iterative 
scheme, we need to compute the boundary condition equa- 
tions and its derivatives with respect to the values of the 
dependent variables at the first meshpoint, the luminosity 
and effective temperature of the star and also with respect 
to the MTR M. 

In view of the need for values of the derivatives with 
respect to the MTR, we have found it very convenient to 
generalize the triangles method presented in Kippenhahn, 
et al. (1967) in the following way. Let us consider a MTR 
and a timestep At; then, the mass of the star will be 
M (1) = M prev + M (1) At. Now, we construct a triangle in 
the HR diagram as proposed by Kippenhahn, et al. (1967), 
i.e., we perform three integrations with 

i) log L/ Lq , log T e ff, , M' 1 ' , 

ii) log L/Lq + A log L/Lq , log T e// , M« , M« , 

iii) log L/Lq , log T e// + A log T eff , M« , M« , 

where logL/Z/Q, logT e // corresponds to one vertex of the 
triangle in the HR diagram and A log L/Lq, AlogT e // arc 
fixed values. We also need to perform integrations with an- 
other MTR value M (2) = M (1) + AM so that, for the same 
timestep, the mass of the star is now M (2) = M prev +M (2) At 
and the triangle in the HR diagram is 

iv) logL/L ,logT e// ,M( 2 >,M< 2 >, 

v) log L/ Lq + A log L/L Q , log T e/ / , M< 2 > , M< 2 > , 

vi) log L/L & , log T eff + A log T e// , , . 

As the model is iterated, M changes. We require the MTR to 



t However, we warn the reader that other authors have employed 
different procedures. See especially, Podsiadlowski, ct al. (2001). 



be in between the values at which the envelopes were com- 
puted, i.e., M (1) < M < M (2) . If not, we change the values 
of M (1) , M (2) . Also, if (log L/Lq, log T e// ), were outside the 
assumed HR triangle we change it with the algorithm pre- 
sented in Kippenhahn, et al. (1967) but now with the same 
mass and MTR values. 

Now, we perform a linear interpolation in M for the 
values of the dependent variables at the bottom of the en- 
velope. Notice that in this way we are automatically inter- 
polating for the correct value of the mass of the star. Also, 
we compute the derivatives of the dependent variables with 
respect to the MTR at each vertex of the triangle. Then, 
all the relevant quantities are found by means of a linear 
bidimensional interpolation inside the triangle in the HR as 
in Kippenhahn, et al. (1967). 

The experience we have gained with our code indicates 
that a very convenient transformation for the MTR variable 



a = In 



M 



IMq/v 



(35) 



Notice that, in this way, the MTR has its correct sign every 
time. The envelopes are computed for given values cti and 
ff2 = cti + Act. For most purposes we have found it enough to 
set Act ~ 0.5, A log L/Lq = 0.02, and AlogT e// = 0.005. 
This represents a good compromise between the amount of 
required envelope integrations and the precision in the in- 
terpolations and allows the whole iteration to converge in 
few steps. 

We compute a set of six outer envelope integration in 
the case that the chemical composition of the outermost 
portion of the star has changed. Otherwise, if the values of 
a, log L/Lq \ogT e ff fall outside the prism corresponding to 
the outer boundary integrations employed for the previous 
model, we change the necessary vertex in order to minimize 
the number of outer integrations. This is important because 
they are time consuming. 

Let us describe how to compute the moment of the on- 
set of mass transfer. We assume that we have begun our 
sequence of models in conditions at which the radius of the 
star is smaller than the corresponding to the sphere equiva- 
lent to the Roche lobe, given by Eq. (25). At each moment 
we compute the structure of the star together with the or- 
bital evolution. Notice that the computation of the orbit is 
necessary in view that there may be angular momentum dis- 
sipation mechanisms operating even without mass transfer 
(e.g. gravitational radiation) . If the components of the pair 
are close enough each other, at some moment the radius 
of the most massive component of the pair-' will overflow 
the Roche radius. At this stage we discard the evolutionary 
model. Then, we guess a MTR low enough and try to con- 
verge the full algorithm presented above. If the MTR value 
guessed is too high and/or the time step is too long we shall 
exchange too much mass. This usually makes the iteration 
divergent, but even if it converges it should be discarded. 
Then, we halve the timestep, and perform a constant mass 



5 Notice, however, that there are situations in which, e.g., a 
canonical 1.4 Mq neutron star has, e.g., a lower mass main se- 
quence companion. In this case, obviously, the fastest evolving 
object will be the low mass one. 
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integration sequence up to Roche lobe overflow. Usually this 
procedure is enough to get a plausible initial mass lossing 
model. From then on, the code works in a very similar fash- 
ion compared to the case of constant mass sequences. This 
procedure for getting the structure of the first model in the 
mass transfer stage is certainly more straightforward than 
the one described by Ziolkowski (1970). 

A more difficult point is to find a self consistent end of 
the mass transfer stage. Usually, we set a low MTR thresh- 
old, much smaller than the value we guessed for the be- 
ginning of the mass transfer stage. As the star gets near 
the end of the mass exchange epoch, the MTR goes down, 
sometimes very suddenly. Meanwhile the MTR is above the 
referred threshold, we assume that the next model will be 
still loosing mass, otherwise we perform a constant mass evo- 
lutionary step and then we compare the radius of the star 
with that of the Roche lobe. If the star is smaller, we assume 
that the mass exchange stage has finished and from there 
on, we continue the evolutionary calculation with constant 
mass. Otherwise, we discard the mass constant integration 
and perform a mass exchanging integration with a halved 
timestep. We have found it to work in most situations. 



6 OTHER RELEVANT NUMERICAL DETAILS 

Another critical, and certainly non trivial problem is the re- 
zoning. This is unavoidably necessary during the course of 
the evolution. Here, rezoning is more problematic compared 
to the case of the evolution of isolated stars. Notice that, the 
stellar radius represents a global characteristic of a star that, 
in the case of isolated stars it appears only in the bound- 
ary condition relating effective temperature with luminosity. 
However, in the case of binary evolution during mass trans- 
fer episodes, it also appears in the condition that equates 
the radius of star with the radius of the Roche lobe. This 
equation, in turn, largely determines the MTR, a quantity 
that deeply affects the structure of the whole star. 

At rezoning we unavoidably introduce numerical noise 
that affects the radius of the star. Thus, it is not surpris- 
ing that, in some critical stages of mass transfer, rezoning 
spoils the convergence of the code. In most of these cases we 
recovered convergence simply performing no rezoning. Even 
in well behaved cases, the noise in the radius of the star 
produces large oscillations in the value of the MTR. This 
is especially true in the case of donor stars with convective 
envelopes. 

In the present version of the code we have simply con- 
sidered the rezoning criteria described in Kippenhahn, et al. 
(1967): We simply test if the variation of the independent 
and dependent variables in a zone is larger than a given 
threshold, if so, we add a new meshpoint dividing the old 
zone in two. At the new meshpoint we compute the value of 
the functions by means of linear interpolations. Conversely, 
if the variation of all variables are below another prefixed 
value, we remove a meshpoint merging two zones in one. As 
it will be shown below, this procedure reveals as adequate 
for the purpose of computing the evolution of low mass CBS 
that give rise to the formation of helium white dwarfs and 
also for the case of high mass CBSs. However, in some cases, 
we have been not able to follow the evolution in very fast 
stages, like a born again phenomenon of a « 0.8 M star 



formed from a 5 Mq object in a 3 days period system. We 
presume such problems to be due to the rezoning. 

Notice that the new meshpoints introduced by linear in- 
terpolations are not solution of the structure equations and 
participate in the calculation in the temporal derivatives of 
pressure, temperature, radius and velocity. We expect to im- 
prove significantly the present rezoning scheme by introduc- 
ing a Hermite interpolation for the values of the dependent 
variables at inserting new meshpoints as suggested by Wa- 
genhuber & Weiss (1994). Hermite interpolation in a zone 
can be performed employing a cubic expression whose coef- 
ficients are computed using the values of the function and 
their derivatives at the edges of the zone, (see Wagenhuber 
& Weiss 1994 for more details). We plan to perform these 
changes in the near future. 



7 APPLICATION OF THE CODE TO SOME 
LOW MASS CLOSE BINARY SYSTEMS 

Before describing the calculations we have performed as a 
test of our new binary evolutionary code, we have to describe 
the physical ingredients we have included in the present ver- 
sion. Because of convenience in the process of preparing the 
code, we have chosen to include some very simplified in- 
gredients, mainly in the nuclear description of the problem. 
We have included the equation of state, OPAL opacities, 
conductive and molecular opacities, and neutrino emission 
processes as in Althaus, Serenelli & Benvenuto (2001). We 
have employed a very simple nuclear reaction network in 
which we have considered equilibrium expressions for pro- 
ton - proton, CNO cycles and helium burning cycles given in 
Clayton (1968). We have followed explicitly the abundances 
of 1 H, 4 He, 12 C, and 1(i O employing the standard Arnett & 
Truran (1969) implicit matrix algorithm. We have neglected 
diffusion processes. 

In order to test our numerical scheme we have chosen to 
compute the evolution of three low mass, CBSs that undergo 
non - conservative, Class A and B mass transfer episodes 
from the main sequence to a highly evolved white dwarf 
configuration. 

In all cases we have assumed a = 0.5 and (3 = 0, i.e. 
that the mass lost by the primary star leaves the system 
(no accretion onto the secondary) and half of the specific 
angular momentum is gone away carried by such material. 
In particular we have considered 

1) M = 2.0 M , q = 1.5, P = 1.5 days 

2) M = 2.0 M , go = 1.5, P = 0.7 days 

3) M = 1.4 M , g = 1.0, P = 1.0 days. 

Here, go = M/M2 is the initial the mass ratio of the binary 
components and Po the initial orbital period. The evolution- 
ary track for the primary component of system 1 is shown in 
Fig. 1 and the conditions at the beginning (points labelled 
with odd numbers) and end (points labelled with even num- 
bers) of mass transfer episodes are included in Table 1. The 
star fills its Roche lobe soon after the end of core hydrogen 
burning. From then on, the star begins to transfer mass on a 
timescale of w 10 7 y (see Fig. 2). The star ends the first mass 
transfer episode with a much lower mass value (see Table 1) 
and begins to evolve bluewards to the white dwarf stage. 
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This evolution is stopped, because the hydrogen envelope 
is still thick enough to ignite in a flash fashion, forcing the 
star to innate on a very short timescale. Then, the star swells 
very suddenly and fills the corresponding Roche lobe again. 
Then, the star suffers from a second mass transfer episode 
depicted in the second panel of Fig. 2. In this second episode 
the total amount of mass transferred is very tiny and occurs 
on a very short timescale. As a consequence the Roche lobe 
size is almost constant and the star is forced to evolve along 
a constant - radius line. Now, after the end of the second 
mass transfer episode, the star has a thinner outer hydrogen 
layer and evolves again bluewards to the white dwarf stage 
but again it suffers from a second hydrogen thermonuclear 
flash. Then, the history is repeated: the star suddenly swells 
up and fills the Roche lobe by third time. Mass transfer 
begins and again a very little amount of mass is transferred 
but now in a much shorter timescale compared to the case of 
the immediately previous mass transfer episode. After end- 
ing this third mass transfer episode, the star evolved again 
bluewards to the white dwarf stage and cools down without 
suffering any other flash. 

As it is well known, hydrogen thermonuclear flashes oc- 
cur very near the stellar surface. As the hydrogen rich outer 
layer is very thin, if we want to compute all these stages 
properly, we do need to consider a very fine zoning for these 
critical layers. As previously mentioned, at flash - induced 
mass transfer episodes, a very tiny fraction of the star is 
lost before the star detaches and starts to evolve bluewards 
again. Then, we do also need a very thin zoning in order 
to account properly for the exact amount of hydrogen re- 
maining in the star, as well as the large gravitational energy 
release when the star undergo large radius excursions. This 
is the kind of physical situation that suggested us the conve- 
nience of designing our binary evolution code in the way we 
have presented in the previous sections, especially regarding 
the treatment of the outer layers and the amount of matter 
we allow to be above the outermost finite - differenced zone. 

We remark that, assuming this fully implicit scheme 
including the MTR we have found that the most difficult 
quantity to compute is just the MTR. This is especially 
true when the outer layers become in convective equilib- 
rium. This is reflected in some strong oscillations shown in 
the first panel of Fig. 2. In any case we feel this not to affect 
the main course of the evolution. Remarkably, such oscilla- 
tions are not present in the flash - induced mass transfer 
episodes. 

System 2 is a Class A system (see Fig. 4) which transfers 
most of its mass when it is still burning hydrogen in the 
core. Then, it suffers three thermonuclear flashes with its 
corresponding mass transfer episodes in a way very similar 
to the corresponding to System I (see Table and Figs. 5 and 
6). 

System 3 (Class B) suffers from four mass transfer 
episodes, three of them induced by each thermonuclear flash. 
The main characteristics of these evolutionary calculations 
are included in Figs. 7-9 and in Tables 2 and 3 we in- 
clude the conditions at the onset and end of mass transfer 
episodes. 

System 3 evolution can be compared with the results 
presented by Podsiadlowski et al. (2001) in their Fig. 14. 
They evolved a system composed initially by a 1.4 Mq nor- 
mal star and a neutron star of the same mass. Mass transfer 



Table 1. Selected evolutionary stages of system 1 



Point 


Age/10 9 y 


M/Mq 


Log(L/L Q ) 


Log(T eff ) 


P[d] 


1 


8.6344152 


2.00000 


1.470297 


3.862974 


1.500 


2 


9.0981136 


0.23161 


0.950006 


3.833146 


2.442 


3 


9.6587255 


0.23161 


2.026359 


4.102230 


2.442 


4 


9.6587278 


0.23096 


2.306213 


4.171956 


2.449 


5 


9.9650518 


0.23096 


1.740276 


4.030472 


2.449 


6 


9.9650530 


0.22935 


2.608005 


4.246795 


2.468 


Table 2. Selected evolutionary stages of system 2 


Point 


Age/W 9 y 


M/Mq 


Log(L/L & ) 


Log(T eff ) 


P[d] 


1 


0.42652401 


2.00000 


1.280174 


3.925827 


0.697 


2 


9.68186661 


0.18421 


0.179932 


3.785182 


1.006 


3 


10.1453763 


0.18421 


1.100802 


4.015836 


1.007 


4 


10.1453778 


0.18395 


1.163948 


4.031070 


1.008 


5 


10.1883408 


0.18421 


1.099444 


4.015355 


1.008 


6 


10.1883427 


0.18348 


1.341075 


4.075097 


1.011 


7 


10.2173678 


0.18348 


1.217356 


4.044674 


1.011 


8 


10.2173692 


0.18313 


1.417843 


4.094100 


1.013 



begins near the end of core hydrogen burning. In doing so, 
they employed a full nuclear reaction network, convective 
overshooting, and different assumptions about mass trans- 
fer from the ones we have done here. In particular, they 
used P — 0.5, i.e., that half of the mass lost by the nor- 
mal star is accreted by the neutron star as long as the rate 
is smaller than the Eddington limit (~ 10~ 8 Mq/?/). Also, 
they considered a = 1 for the angular momentum losses. 
Clearly, these are large differences compared with our Sys- 
tem 3 and also with the physical ingredients we considered in 
the present paper. Notice that Podsiadlowski et al. (2001) 
found three hydrogen thermonuclear flashes, two of which 
leading to Roche lobe overflow. In our system, we also have 
three hydrogen flashes, all of them leading to lobe overflow. 
This is an important difference that we suspect it to be re- 
lated with the different assumptions in the orbital evolution 
of the system. We also consider that the differences in the 
treatment of nuclear energy release are of key importance. 
Clearly, our equilibrium cycles predict stronger flashes. How- 
ever, notice that the mass of the final helium white dwarf is 
very similar in both calculation (0.199 Mq for Podsiadlowski 
et al. 2001). 

Finally, in Fig. 10 we show the evolution of the central 
part of the three computed objects in the temperature - 



Table 3. Selected evolutionary stages of system 3 



Point 


Age/ltfy 


M/Mq 


Log{L/L Q ) 


Log(T eff ) 


P[d] 


1 


2.0641016 


1.40000 


0.844531 


3.797418 


1.000 


2 


3.4962904 


0.20567 


0.472654 


3.702205 


2.803 


3 


3.8905910 


0.20567 


1.053372 


3.847385 


2.803 


4 


3.8905929 


0.20530 


1.278767 


3.903560 


2.809 


5 


3.9239877 


0.20530 


1.299815 


3.908821 


2.809 


6 


3.9239887 


0.20490 


1.546744 


3.970370 


2.815 


7 


4.0033926 


0.20490 


1.385939 


3.930169 


2.815 


8 


4.0033931 


0.20415 


1.853513 


4.046706 


2.827 
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density plane. As expected, the evolution is rather similar 
and is limited to temperatures below helium ignition. The 
flash episodes suffered by these objects are responsible for 
the oscillation in the temperature of the star before it gets 
near the final density value. Then, the objects cool down 
almost at constant central density, as it should be expected. 

These calculations indicate that the above described 
scheme is capable of handling the evolution of close binary 
systems in a very adequate way. In our opinion, the here 
described code should be a valuable tool in handling binary 
evolution. Let us remark that, we have tested the code with 
low mass systems because are the ones we are most inter- 
ested in. Also, we have found very good convergence of the 
method for the case of massive CBSs. 

In closing, we should remark that the above described 
results should not be considered a state - of - the - art evo- 
lutionary calculations but a serious test of the numerical 
scheme as a previous step for future works. This is so mainly 
because of the lack of a detailed nuclear reactions network 
and also for the absence of a treatment of diffusion. Notice 
that diffusion is a key process in the evolution of objects 
in the range of masses we have here considered and even 
largely determine the total number of thermonuclear flashes 
the object suffers from (see Althaus, et al. 2001). 



8 DISCUSSION AND FUTURE 
APPLICATIONS 

In this paper we have presented a Henyey-type code tai- 
lored for computing stellar evolution in close binary systems 
(CBSs). This code is a modification of the scheme presented 
long ago by Kippenhahn, Weigert & Hofmeister (1967) for 
the case of the evolution of isolated objects. 

Here we considered that mass transfer occurs when the 
star overflows its Roche lobe. The main characteristic of the 
present code is that it automatically computes the beginning 
and end of mass transfer stages, computing the MTR of the 
object in a fully implicit, self consistent fashion. Perhaps, 
main shortcoming of the scheme we have presented is that 
it is not able to compute common envelope evolution, simply 
because at these stages it is not possible to equate the radius 
of the star to that of the Roche lobe. 

In the present version of the code, we have considered 
radiative OPAL and conductive opacities, neutrino emissiv- 
ities, convection following the mixing length prescription, 
and a detailed equation of state. The nuclear reactions have 
been considered only by means of the equilibrium cycles. 
Regarding the physics of mass transfer, we have considered 
the possibility that a fraction of the mass transferred from 
the primary is accreted by the secondary. Also we have in- 
cluded angular momentum dissipation due to gravitational 
radiation and magnetic braking. 

We have found that the scheme detailed above is able to 
handle most of the evolutionary stages suffered by a donor 
star belonging to a CBS. Apart from the results we have 
presented as a first application, we have also tested it for 
the case of massive CBSs finding it to be very adequate also 
for these purposes. Generally speaking, the code is able to 
perform single stellar evolution and also to find the onset 
and end of each mass transfer episode. In all the test cal- 
culations we have performed we found it much more easy 



to compute a mass transfer rates (MTRs) when outer layers 
are in radiative equilibrium. In the case of convective lay- 
ers convergence is attained but sometimes it is necessary a 
larger number of iterations. 

As a first application of this code to astrophysically rel- 
evant conditions, we have applied it to the case of the forma- 
tion of low mass helium white dwarfs. These computations 
should not be considered as a state - of - the - art ones 
but as a test of the numerical scheme as a previous step for 
future works. This is so mainly because of the lack of a de- 
tailed nuclear reactions network and also for the absence of 
a treatment of diffusion. 

We expect that the physical and numerical performance 
of the code should be largely improved if we replace the outer 
boundary condition for the stellar radius by the formulation 
of Ritter (1988) for mass loss transfer (see the Introduc- 
tion). This should be especially important at the beginning 
and end of each mass transfer episode and/or in the case of 
very low mass transfer rates. Also, the implementation of a 
moving grid should help in producing smooth curves of mass 
loss vs. time. 

As a future application for the present code we plan to 
study the problem of the formation of helium white dwarfs 
including the complete physical ingredients we have consid- 
ered in the previous works of our group on this topic. Also, 
we plan to study the relevance of the irradiation (Tout, et al. 
1989; D'Antona & Ergma 1993) in the evolution of the mem- 
bers of CBSs. This should be especialy relevant for the case 
of LMXBs which in some evolutionary stages are expected to 
undergo mass transfer driven by irradiation (Podsiadlowski 
1991). 

Fig. 4-10 are available upon request to the authors at 
their e-mail addresses. 
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Figure 1. The evolutionary track for the primary component of system 1 (M = 2.0 Mq , go = 1-5, Po = 1.5 days). Dots labeled with odd 
(even) numbers indicate the onset (end) of mass transfer episodes. Notice that the star undergoes two hydrogen thermonuclear flashes 
each of them producing the large loops shown. For more details, see text. 
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Figure 2. The mass transfer rate from the primary star of system 1. The star undergoes a long mass transfer episode soon after core 
hydrogen exhaustion (panel A). In panels B and C we depict the mass transfer rate corresponding to the episodes due to hydrogen 
thermonuclear flashes. Notice that they are, by far, shorter and with a maximum rate about two to three orders of magnitude larger 
than the corresponding to the first episode. In each panel, time is counted from the onset of mass transfer on. For the age of the star at 
these moments, see Table 1. 
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Figure 3. The Orbital period and semiaxis as a function of the mass of the primary star of system 1. 
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